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ABSTRACT 



The nonequilibrium dynamics of the chiral phase transition 
expected during the expansion of the quark-qluon plasma pro- 
duced in a high energy hadron or heavy ion collision is studied 
in the 0(4) linear sigma model to leading order in a large N ex- 
pansion. Starting from an approximate equilibrium configuration 
at an initial proper time r in the disordered phase we study the 
transition to the ordered broken symmetry phase as the system 
expands and cools. We give results for the proper time evolu- 
tion of the effective pion mass, the order parameter < o > as 
well as for the pion two point correlation function expressed in 
terms of a time dependent phase space number density and pair 
correlation density. We determine the phase space of initial con- 
ditions that lead to instabilities (exponentially growing long wave 
length modes) as the system evolves in time. These instabilities 
are what eventually lead to disoriented chiral condensates. In 
our simulations, we found that instabilities that are formed dur- 
ing the initial phases of the expansion exist for proper times that 



are at most 3/m/c and lead to condensate regions that do not 
contain large numbers of particles. The damping of instabilities 
is a consequence of strong coupling. 
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1 Introduction 



In recent years there have been numerous investigations concerning the possi- 
bility of forming large correlated regions of misaligned vacuum during highly 
energetic collisions |], @, |3|. Those regions, in which the quark condensate 
< qi q.j > is nonzero but points along the wrong direction in isospin space 
have been named Disoriented Chiral Condensates (DCC's). If large DCC's 
are indeed formed, they would produce spectacular events in which one could 
observe strong correlations between the emitted pions. In fact, the original 
idea of DCC's was invented [0 to explain rare events observed in cosmic 
ray experiments M: the Centauro events, where there is a deficit of neutral 
pions. In this approach the explanation of the deficit would be the result of 
the decay of domains in which the condensate has a vanishing component in 
the 7T° direction. Anti-Centauro events (whose experimental status is still 
uncertain), where the emission would be predominantly neutral would be 
explained as coming from regions in which the condensate points along the 
7T° direction. 

If these events are a result of the creation of DCC's, this would be di- 
rect evidence for the existence of a chiral phase transition in the plasma 
formed following an ultrarelativistic collision and would allow us to explore 
the physics of the chiral phase transition. The possibility of producing DCC's 
in high energy collisions has originated several experimental proposals P, f| 
and a number of theoretical papers - fL2| . It has been recognized that 
the possibility of forming large regions of DCC relies on the existence of a 
substantially large regime in which the hot plasma formed after the collision 
evolves out of equilibrium ||. In fact, if thermal equilibrium is approxi- 
mately preserved by the dynamics, the typical correlation length would be 
determined by the pion mass and therefore would be too small to matter. 
For this reason, there have been a number of authors making different at- 
tempts to analyze the nonequilibrium aspects of the dynamics of the chiral 
phase transition. These attempts vary in form and content: some authors 
performed numerical simulations on classical models || ||, others used 
phenomenological terms - inspired in classical kinetic theory - to model the 
interaction between the condensate and the quasiparticles while some 



attempted to incorporate quantum and thermal fluctuations || [12] into the 
theoretical framework. 

There are clearly two important questions that any theoretical model 
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should answer. The first one is to determine whether during the evolution 
that follows the collision there are instabilities affecting the fluctuations. If 
this happens,then there is a chance for the correlations to grow. If this occurs, 
structure may tend to form through a process like spinodal decomposition. 
In our simulations of the chiral transition in the sigma model we find that 
several types of initial fluctuations involving the time derivative of the a field 
lead to instabilities, so that we indeed find that the first ingredient is present 
in our simulations. 

The second question is if, assuming the instability exists, the correlated 
domains can grow large enough so that many pions can be emitted from 
each domain making the detection of DCC's possible. The time scale for the 
instabilities to exist is strongly influenced by the strength of the interactions. 
In fact, some of the intuition one may have developed by analyzing similar 
problems in other contexts, such as the cosmological one, where the coupling 
is extremely weak, may not apply at all. Using the linear sigma model as the 
effective field theory valid near the chiral phase transition temperature, we 
find that we need to be in a regime of strong coupling (A r ~ 10) in order to fit 
low energy scattering data. In the strong coupling regime instabilities exist 
only for a short period of proper time as we shall see below. This makes it 
difficult to obtain large sizes for the correlated domains. On the other hand 
in a weakly coupled system the domains can grow for long periods of time 
and a totally different picture can emerge. 

In our approach we do not put in instabilities by hand but look at typical 
fluctuations of initial conditions starting in a region of stability. Previous 
researchers have assumed the existence of an initial instability and focused 
their calculations on studying the growth of correlations. The typical and 
quite drastic way in which the instability had been previously introduced in 
the problem was by using the so called quench approximation. In a quench, 
the system (the hot plasma created after the collision) is subject to a sud- 
den external action (the expansion into the vacuum) that has the following 
net effect: it does not change the state of the system, but instantaneously 
modifies the effective potential turning it "upside down". This is clearly an 
idealization that may be entirely inappropriate. In this paper we will present 
an approach that does not require this approximation and that enables us 
to study (in a simplified model) the existence of instabilities in a rather 
straightforward way. 

Thus, our aim is to study the evolution of DCC's without imposing the 
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quench approximation or any other ad hoc way of modelling the appearance 
of an instability. That is we start in a situation where the particle masses 
are positive and one is in the symmetric phase. We start at an initial value 
of the proper time Tq in the symmetric phase by choosing a thermal distri- 
bution of particles above the critical temperature. This initial condition is 
not necessary, but is one way of ensuring that the initial system is in the 
disordered phase. We imagine that the system cools via the expansion of 
the system into the vacuum. By cooling we mean a reduction in the energy 
density with proper time. If we were in equilibrium, the temperature would 
automatically decrease with the proper time. The simplified picture we have 
of cooling is the boost invariant picture popularized by Bjorken fl3fl . This 
picture is consistent with various hydrodynamic approaches to hadronic as 
well as heavy ion collisions. We kinematically constrain expectation values 
to only depend on the fluid proper time r = y/t 2 — x 2 , where x is the distance 
along the collision axis in the center of mass frame. The natural expansion 
and subsequent cooling of the plasma into the ordered vacuum is studied 
numerically by solving the update equations for the quantum modes as well 
as for the proper time evolving expectation values. In this way, we are able 
to study the evolution of the plasma in a self consistent way without impos- 
ing any instability by hand. We analyze various reasonable initial conditions 
on the fields and determine whether they lead to instabilities. That is for 
various initial conditions we determine the proper time evolution of the effec- 
tive pion mass. When the effective pion mass becomes negative instabilities 
ensue. We also determine the time evolving order parameter < a > and 
the adiabatic phase space number density and pair density which determine 
the spatial correlation function for the pion field. These number densities 
are related to physical measurables (such as the rapidity distribution of fi- 
nal particles) at later times. The initial conditions which lead to the largest 
instabilities have initial velocities in either the a or ir directions. We com- 
pare the nonequilibrium results for the number density with those that would 
have resulted from an expansion with local thermal equilibrium in the co- 
moving frame. When instabilities arise the distributions tend to narrow in 
momentum space, especially in the transverse direction. 

In the investigations done so far, simple phenomenological models have 
been used hoping that they describe the fundamental physics involved in the 
dynamics of the phase transition. We will employ the linear sigma model, 
the most popular one in this context, which seems to have the essential 
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attributes of being simple but realistic enough: it appropriately describes the 
low energy phenomenology of pions and has also the correct chiral symmetry 
properties. The initial conditions we will impose are motivated by matching 
it to the situation one expects to attain in a highly energetic collision. We will 
assume that the quantum state of the system is a thermal density matrix at 
an initial instant of proper time Tq. We will choose the initial temperature T 
to be slightly above the critical temperature for being in the disordered phase. 
As shown in the appendix the critical temperature is given by: T c 2 = 3/ 2 . We 
will choose the parameters of the model to give reasonable values for three 
experimentally determined quantities : the mass of the pion m n , the pion 
decay constant f n and the s wave n — n phase shifts above threshold. These 
three measurements completely determine the parameters of the model. One 
important constraint on this model is the triviality of the model as the cutoff 
is removed. The theory only makes sense at cutoffs below the Landau pole 
which occurs at a value of the cutoff A when the bare coupling constant first 
becomes negative for positive renormalized coupling constant. This limits 
the size of the renormalized coupling constant. The maximum renormalized 
coupling constant as a function of A obeys for large A 
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Since the mass difference between the a and 7r is directly proportional to X r , 
this leads to an upper bound for the the mass of the a resonance as a function 
of A. Therefore, unlike at tree level, the mass of the a in the fully quantized 
theory is constrained in this model. A is also constrained from the physical 
considerations that we want the mass of the o to be less than the cutoff. 
However the mass of the a resonance increases as we decrease the cutoff 
since then the renormalized coupling increases. This pins down the cutoff 
to lie between 700 MeV and lGeV. The rest of the paper is organized as 
follows. In section 2 we describe the linear sigma model in the leading order 
in the 1/N expansion. In section 3 we discuss the Baked Alaska Scenario 
and we derive renormalized update equations for the proper time evolution 
of the field theory. In section 4 we describe the initial conditions we use 
to study the development of disoriented chiral condensates. In section 5 we 
discuss the results of extensive numerical simulations. In the appendix we 
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rederive for completeness the properties of the linear sigma model in the large 
N approximation. 



2 Model and approximations 

We will use the 0(4) linear sigma model to describe the evolution of the 
pions. We are well aware of the limitations of this approach, that provides 
a reasonable phenomenological model only for a limited range of energies 
(typically smaller than 1 GeV). In spite of its shortcomings, this model 
captures some of the essential physics involved in the dynamics of the phase 
transition that may produce disoriented chiral condensates. In particular the 
chiral phase transition takes place at a reasonable temperature (T c = \Z3f w ), 
and the low energy n - n scattering amplitudes are reasonable in this model. 
The mesons are organized in an 0(4) vector $ = (a, fr) and the action is (in 
natural units h = c = 1) 

S = J d A x{)^d<S> • <9$ - -A($ ■ $ - v 2 ) 2 + Ha}. (2.1) 

where we have use the Bjorken and Drell metric: (1,-1,— 1,-1). We will 
describe the evolution of the mean value $ =< $ > and the two point correla- 
tion functions including the effects of quantum and thermal fluctuations. For 
this purpose, we will treat ( |2.1| ) as a quantum theory defining initial values 
(i.e., the quantum state of the system at a given proper time) and will follow 
the real proper time evolution of $ and the correlations. The complexity of 
the problem forces us to adopt some approximations. Perturbation theory 
is useless for our purpose [[14], [l5j and a scheme which is non-perturbative 
in A must be adopted. In this context, the most popular approaches which 
allow for a real time analysis are the the Hartree (or Gaussian) ansatz |L6] 



and the large N expansion of the 0(N) sigma model [M. We will adopt 
the later since it presents several advantages. On the one hand, the expan- 
sion is systematic and allows us to study higher order corrections (work is in 
progress in this direction and results will be presented elsewhere (ll|). On 
the other hand, when using the Hartree ansatz in the context of the study 
of DCC's one is forced to take also the large N limit (see Ref. fl2]l). This 
is due to the well known fact that the Gaussian approximation violates the 
Goldstone theorem giving an unphysical (and not necessarily small) mass to 



6 



the pions in the H = limit. Of course, the approximation we adopt here 
is not expected to capture all the features of the phase transition (as is well 
known, mean field theory fails to predict the correct critical exponents but 
allows us to explore the strong coupling regime). Because of the triviality of 
the O(N) sigma model as one takes away the cutoff, an aspect of the exact 
theory that is preserved in the large N approximation, we have to seriously 
take into account the cutoff and its ramifications. One of these ramifications 
is that the renormalized coupling constant has a modest upper bound of the 
order of ten at a cutoff of one GeV. This gives an upper bound to the mass 
of the a resonance whose value depends on our choice of cutoff. 

The large N effective equations can be obtained in a variety of ways, which 
are extensively discussed in the literature [171]. A very convenient method is to 



use an effective action, which is a functional of the mean values of the original 
fields $ and of an auxiliary constrained field x 0- We start with a classical 



action S"[$,x] constructed from ( |2.1| ) by replacing x — A($ 2 /2iV — v 2 ). As 
this action is now quadratic in $ we can perform the functional integral 
over those fields and are left with a functional integration over x which, to 
leading order in 1/N can be calculated by the stationary phase method. In 
the appendix we review the details of this calculation and calculate all the 
propagators and the n — n scattering amplitude in the leading order in the 
large N expansion. Higher order corrections can be systematically computed 



in this way [19[ and an expansion of the effective action r[$, x] m powers of 



1/N can be obtained [ 17fl . We will consider here only the leading order terms 
which give the following equations (for notational convenience we drop the 
overbars and denote the expectation values $ and x): 

(n x + x {x))$i{x) = H5 tl (2.2) 
x (x) = \(-v 2 + $ 2 (x)+NG (x,x)) (2.3) 
G \x,y) = i{n x + x (x))S{x-y). (2.4) 

The structure of these equations is indeed very simple. The field x plays 
the role of the effective mass for the mean values <I>j and satisfies the "gap 
equation " Q2.3|) . The function Gq(x,x) that appears in fl2.3|) is the coinci- 



dence limit of the propagator G (x, y) that inverts the operator G 1 defined 
in Q2.4). We can use an auxiliary quantum field <p(x) where < <p(x) >= 



and 

(n x + x(x)) ( f ) (x) = (2.5) 
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to determine Go(x, y).We construct the propagator as 

G {x,y) =< T (f)(x)(f)(y) > 

. 5ijG (x,y) is the the pion propagator when < ir^x) >= 0. The initial 
value problem associated with equations ([2T2|)-([2T4|) will be solved in the 
next section. Here, we would like to address the issue of how to use this 
model to make contact with the phenomenology we want to describe. Thus, 
we must fix the values of the parameters appearing in the above equations 
so that they describe low energy pion physics. The measurable quantities we 
want to reproduce are the pion mass = 135 MeV, the pion decay constant 
f n = 92.5 MeV as well as the s wave, 1 = phase shifts in the energy range 
300 — 420MeV. To fit these physical quantities we analyze our equations in 
the "true vacuum" state (i.e., in equilibrium at zero temperature). In such 
state the derivatives of the expectation values vanish and we have $ = (o~„, o), 
X — Xv where a v and Xv are some constants whose values will be determined 
below. The physical masses can be related to the parameters of the theory by 
computing the inverse propagators of the pion and sigma fields. The s-wave 
phase shifts is determined from the it — it scattering amplitude obtained in 
this approximation which is given by the exchange of the composite field x 
propagator in all three channels. A complete review of the vacuum and finite 
temperature properties of this model in leading order in large N is found in 
the appendix and we summarize the results here. The vacuum expectation 
value of a is determined by f n 

<J V = U- (2.6) 

On the other hand, in the vacuum we have 7? = and the pion inverse 
propagator is Gz\ n (x, y) = Gq 1 (x, y)Sij. Therefore, the vacuum expectation 
value of x is 

Xv = ml = m 2 . (2.7) 

The sigma mass can be approximately determined in terms of the inverse 
sigma propagator as the zero in the real part of the inverse propagator (a 
more precise determination which gives a slightly different result is from the 
peak in the 1 = 0,1 = scattering amplitude). This leads to the equation: 

m 2 a = m 2 + a 2 Re[G 0xx (ml)}, (2.8) 
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where Go xx ^ s the composite field propagator in the absence of symmetry 
breaking, i.e. 



6 o> 2 ) = ^ + yn( P 2 ), 



(2.9) 



and the polarization II = iG^ is given by 



n( P 2 ) = -iJ [A] (x - q 2 r\x - (p + q?r\ 



(2.10) 



and is explicitly given in the appendix. We use the notation that [d n q] = 
d n q/(2n) n . In the absence of symmetry breaking the composite field propa- 
gator is the geometric sum of the bubble chains. In the presence of symmetry 
breakdown one is also summing the graphs where the bubble is replaced by a 
single propagator and two tadpoles (this includes the sigma exchange graphs). 
One can determine the bare mass of the pion in terms of the physical pion 
mass using the gap equation, since Xv = m 2 



This relationship tells us that the bare mass goes to negative infinity 
quadratically with the cutoff, as is true for the exact lattice theory. (That is 
the reason why in strong coupling the theory is related to the Ising Model). 
This relationship also allows us to eliminate the bare mass from the theory 
as well as the quadratic dependence on the cutoff. The logarithmic depen- 
dence on the cutoff are removed by coupling constant renormalization. We 
obtain a value of the bare coupling A at a fixed cutoff A (or equivalently the 
renormalized coupling A r ) by comparing our large-N result with the Pade fit 
of Basdevant and Lee |22] to the 1 = s-wave scattering amplitude. This is 
discussed in detail in the appendix. Once A is determined both v and m a are 
fixed. Therefore all the bare parameters of this theory are fitted in terms of 
the three experiments which determine m, f w , and A. The external field H is 
determined from the equilibrium time independent solution to field equation 



m 2 = -\ v 2 + \fl + \NG (x, x) 



(2.11) 



where 




(2.12) 
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for the a field: x a = m nf^ = H. Our numerical simulations automatically 
have a cutoff since we are performing numerical simulations in a box. Thus 
one can use either the bare or renormalized parameters to describe the prob- 
lem. The renormalized parameters, however are determined from physical 
measurements and are more fundamental. However since the a model is not 
an asymptotically free field theory, it is a trivial theory when we take the 
cutoff away. Thus one is not allowed to take the continuum limit and must 
consider the theory as an effective field theory which is not valid for energies 
in excess of a few GeVs. Our model is a theory which makes physical sense 
(such as having a well defined ground state) only with a physical cutoff A p 
(see |20| for discussions on $ 4 as a cutoff theory). On physical grounds 



we want 2m n < m a < A. This relation is quite a constraint on A since when 
we lower A we increase the value of the mass of the sigma. We do not expect 
the theory to be valid for energies above 1 GeV since in that regime the 
correct dynamics is described by QCD. In fact if we raise A above 1 GeV, 
the maximum value of the renormalized coupling goes down and the a mass 
becomes unacceptably low. Reasonable values for the mass of the a constrain 
A to be in the range .7GeV < A < lGeV. In that range the best value of 
the renormalized coupling A r is between 7-10. For this range of values for 
A r this model agrees qualitatively with low energy scattering data. With 
this choice of parameters, it is difficult to obtain values of the a mass higher 
than 450MeV. Based on the perturbative calculations of Basdevant and Lee 



22j which were then subjected to a Pade improvement, and the connection 
between 1/N expansions and resummations we expect that at next order in 
1/N this upper bound on the a mass will be raised slightly. 

As will be clear from our discussion below, when solving the equations we 
can use a different value of the cutoff provided we scale the bare couplings 
appropriately (so that we keep the physical quantities unchanged). As the 
theory we are using is renormalizable, this scaling is well defined and known. 
However, the cutoff can not be taken too large since the theory only makes 
sense for positive values of the bare coupling. This is the Landau pole prob- 



lem and is related to the triviality of the theory. |2Bj ]23| . As is discussed in 
detail in the appendix, the bare and renormalized couplings are related by 
the equation: 

A = — — t— r. (2.13) 

1 - NA r n(0) v ; 
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We see that A has a pole when 1 = iVA r II(0). At large values of the cutoff 
this Landau pole occurs when 

log(2A/m 7r ) » (2.14) 

For fixed A r we notice that for values of the cutoff larger than that given by 
( [2.141) ) the bare coupling becomes negative which makes the theory undefin- 



able as a Euclidean lattice field theory. This relationship also means that 
at a fixed cutoff, there is a maximum A r that we can consider. That is the 
inverse relationship is 

Ar = i + Nxmo) ■ (2 - 15) 

The maximum of A r is reached at infinite A and one obtains: 

x ™=Nwy (2 - 16) 

For a cutoff of lGeV the maximum A r is 13. 

As a final comment we would like to point out that the effective action 
method we employed can also be used to compute higher order corrections in 
1/N. In that case it is necessary to use a formalism that enable us to derive 
real and causal equations for the expectation values (otherwise one cannot 
even pose the initial value problem). This formalism is known as the "closed 
time path" method developed by Schwinger, Mahanthappa, and Keldysh JZ5 



The analysis of the 1/N corrections obtained using this approach (that carry 
the effects of the collisions that may produce relaxation towards equilibrium) 
will be analyzed elsewhere |1| (see also Ref. |6|, [14], [lTj for recent applications 
of this formalism in the DCC and other related contexts). 



3 Cooling by expansion and Baked Alaska 
scenario 

3.1 The basic idea 

To analyze the possibility of forming DCC's we should take into account the 
specific features characterizing the situation after a highly energetic heavy 
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ion collision. Experimentally, a flat plateau in the distribution of produced 
particles per unit rapidity is observed in the central rapidity region (these 
results are obtained in p-p and other collisions). This suggest the existence 
of an approximate Lorentz boost invariance. Thus, the simplest picture of 



a collision, due to Landau p6[, is one in which the excited nuclei are highly 
contracted pancakes receding away from the collision point at approximately 
the speed of light. The boost invariance implies that the evolution of the 
"hot plasma" that is left in between the nuclei looks the same when viewed 
from different inertial frames. Of course, this is an approximate picture that 
is not valid for large values of the spatial rapidity and for transverse distances 
of the order of the nucleus size. The existence of this approximate symmetry 



can be used to make a very simple hydrodynamical model [|I3] that, in some 
cases, may describe the evolution of the plasma. It is worth reviewing very 
briefly these ideas. First one should recognize that the natural coordinates to 
make a boost invariant model are the proper time r and the spatial rapidity 
rj defined as 

r^{f-x^\ ^llog(^). (3.1) 

The observed symmetry will be respected by the model if one imposes initial 
values on a r =constant hypersurface (and not at constant laboratory time 
t). 

If we had local thermodynamic equilibrium during the expansion so that 
the relaxation rate is faster than the expansion rate and assumed that we 
could approximate the field theory dynamics with a hydrodynamic flow, one 
would find by solving the hydrodynamical equations with this type of initial 
conditions (homogeneity along the constant r surface) that the energy density 
drops as r~ a . Here a = 1 + Cq where Co is the speed of sound, which depends 
on the equation of state of the fluid p = c^e. In the ultrarelativistic case 
Cq = 1/3 and the temperature falls as r^ 1 / 3 . 

According to this simple model of a collision, the plasma evolves in a 
highly inhomogeneous way when viewed from the laboratory frame. In fact, 
analyzing a constant t surface we realize that the field configurations strongly 
depend on the spatial coordinate x. Near the light cone |x| = t the system is 
"hot" (corresponding to small values of the proper time r). On the contrary, 
for small values of x (that correspond to larger values of r) the system is 
"colder" . This type of configuration, hot in the outside and cold in the inside, 
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is schematically known as Baked Alaska (or a Swiss Bombon, according to 
other culinary traditions). 

In this paper we want to study the formation of DCC's using some of 
the ideas presented above. We will not assume a quasi equilibrium situation 
or use a phenomenological hydrodynamical (or kinetic) model to describe 
the evolution of our system. On the contrary, we will study the nonequilib- 



rium evolution in its full glory and solve equations (|2.2| )- (|2.4j ), which include 
thermal and quantum effects. Using the coordinates (|3.1|) and fixing boost 
invariant initial conditions at an initial proper time tq (whose numerical value 
we discuss below) we introduce the expansion and "cooling" of the plasma 
in a natural way. Therefore, we do not need to introduce any ad hoc cooling 
mechanism by hand. The cooling, if any, will appear as a result of the evo- 
lution, which is fully out of equilibrium. In this way, we can really test the 
validity of the "quenching" approximation that has been almost universally 
used when analyzing the evolution of DCC's. 

3.2 The equations 

Our treatment is very similar to the one required to study quantum field 



theory in a curved spacetime [^J . Thus, in the coordinates ( |3.1[ ) Minkowsky's 
arc element is 

ds 2 = dr 2 - T 2 drf - dx\, (3.2) 

which has the same form of a Kasner Universe (an anisotropically expanding 
Universe, which in this case is nevertheless flat, see ||27|| ). To study the 



evolution of the mean fields and correlations in this coordinates our first task 
is to rewrite equations ( |2.2| )-( ]2"~4] ) using the new variables. Assuming that 



the mean values $ and x are functions of r only (homogeneity in the constant 
t hypersurface) we have 

r- l 8 T rd T Mr) + x(r) ^(r) = HS a (3.3) 

x ( r ) = \(-v 2 + $, 2 (r) + N < <f> 2 (x, r) >) (3.4) 
where the quantum field 0(x, r) satisfies the Klein Gordon equation 

(r" 1 ^ rd T - r- 2 d 2 - d\ + x{x))<f>{x, r) = 0. (3.5) 
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The quantum field <j>(x, r) denned here is an auxiliary field which allows us to 
calculate the Wightman function Go(x,y) by taking the expectation value: 

G (x, y; t) =< cj)(x, r) <j>(y, r) > . (3.6) 

When the expectation value of the pion field is zero, then this field corre- 
sponds to one component of the pion field. 

As usual, we expand this field in an orthonormal basis 

<f>{r),x±, r) = J[d 3 k](exp(ik • x)/ k (r) a k + h.c.) (3.7) 

where k ■ x = k v r] + k±x±, [d 3 k] = dk v d 2 k±/ (2n) 3 and the mode functions 
/k( r ) evolve according to (a dot here denotes the derivative with respect to 
the proper time r): 

A + (^ + £ + x(t) + ^)A = 0. (3.8) 

For the creation and annihilation operators to obey the usual canonical com- 
mutation relations [cife,a|] = 1, then the mode functions must satisfy the 
Wronskian condition: 

//• - /•/ < 

The expectation value < (f> 2 (x, r) > can be expressed in terms of the mode 
functions / k and of the distribution functions 

n k =< a[a k >, g k =< a k a k >, (3.9) 

which entirely characterize the initial state of the quantum field. For sim- 
plicity, we will assume that the initial state is described by a density matrix 
which is diagonal in the number basis (like a thermal state). In such case, the 
only nonvanishing distribution is n k . Thus, replacing the above expressions 
in Q3.4D we have: 

X(r) = X(-v 2 + *?(r) + \n J [rf 3 k]|/ k (r)| 2 (1 + 2 n k )). (3.10) 

We assume the initial density matrix is one of local thermal equilibrium 
in the comoving frame. In the comoving frame, the expectation value of the 
energy momentum tensor is diagonal and of the form 

diag(e,p,p,p). (3.11) 
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We then have at r = r (the surface of constant energy density and temper- 
ature T ) that: 

Hk = e Jl _ \ ' ( 3 - 12 ) 

where A, = 1/T , El = yjk>(r ) + X (t ) and k 2 = k 2 /r 2 + k\. 

We can also use a time dependent set of creation and annihilation opera- 
tors to describe the quantum field 0. If we use the first order adiabatic mode 
function: 

e -iVk{r) 

fk = v /2^ ' ; dyk/dr = uj k) (3.13) 

where Ct^(r) = ( + &i + x( r )) 1 ^ 2 Then if we expand the field in terms 
of these mode functions we h 

0(77,^,r) = ^ j[d 3 k](exp(*kx)/ k °(r) a k (r) + h.c). (3.14) 

We can now define the first order adiabatic interpolating number and pair 
distribution functions via: 

nk(r) =< 4(r)a k (r) >, g k (r) =< a k (T)a k (r) > (3.15) 

The time dependent creation and annihilation operators satisfy 

d/° + d t /*° = (3.16) 

in order that the usual canonical commutation relations hold. Then we can 
rewrite the Green's function G° in the two bases: 

G°(x,y;r) = - f[d\]e^^-\- (1 + 2 n k (r) + 2Re[g k (r)e-^]) 

(3.17) 



or 



G\x,y-T)= l -J[d^y k ^ [l + 2n k (0)] |/ fe (r)| 2 . (3.18) 
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The t dependent variables n and g now have the physical interpretation 
as the interpolating phase space number and pair density for the case when < 
7Tj >= 0. These operators agree with the time independent number and pair 
densities defined by ( |3.9| ) at r , and in the out regime become the physically 
measurable number and pair densities. For < 7Tj >^ the actual pion two 
point function is more complicated than G° and this is only one piece of 
that Green's function. Of course in the vacuum Isospin conservation requires 
< iii >= 0, so that once we are approaching the final true vacuum state 
during the cooling process, one can use these interpolating number and pair 
operators to describe the physics of the problem. Also note that during the 
instability phase, io\ can become negative and for those modes the adiabatic 
basis does not exist because the Wronskian condition for the mode functions 
is no longer satisfied. 

The two sets of creation and annihilation operators are connected by a 
Bogoliubov transformation: 

a k (r) = a(k, r)a k + (3(k, r)a_ k . (3.19) 

a and /3 can be determined from the exact time evolving mode functions via: 

*(k,r)=i(f^-^f k ) (3.20) 

mr)=i(f° k ^-^f k ) (3.21) 

In terms of the initial distribution of particles n k (ro) = n k and the Bo- 
goliubov coefficients a and (3 we have: 

n k (r)=n k + \/3(k,T)\ 2 (l + 2n k ) (3.22) 

g k (r) = a(k,T)P(k,T)(l + 2n k ) (3.23) 
3.3 Initial conditions and renormalization 



We desire to solve equations ( |3.3|) ,(|3~lf) and fl3.10|) as an initial value problem 



To do this we need to give Cauchy data (the function and its derivatives) for 
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the mean values $i(r) and for the mode functions fk{^)- This, together with 
the distribution function and at r fully defines the problem. Notice 
that the initial value of xi T ) is determined self-consistently by solving the 
gap equation ( |3.10| ) at the initial time r . 



The quantum state is fully determined by the distribution and by the 
initial data f^(ro), /k(To)- These data fix the vacuum state upon which the 
Fock space is built. The definition of the vacuum state in a curved spacetime 
is certainly a tricky issue with a rather long history p7 |. Fortunately this is 



not an important problem for us, due to the fact that the background space- 
time is flat. Thus, it is possible to choose the initial data fk(T ) and /k(To) 
so that the vacuum state coincides with the ordinary Minkowsky vacuum, 
at least for high momentum. It can be shown that this is accomplished by 
taking the "zero order adiabatic" vacuum where 

where lo^t) = (k 2 /r 2 + k\ + x(t)) 1 / 2 . We will fix these initial conditions 
for the normal modes, which are normalized according to the Wronskian 
condition /*/ — /*/ = i. The distribution will be taken as thermal, 
characterized by a temperature T: 

n k = [eM^(r )/k B T) - l]" 1 . (3.25) 

As we said before, the existence of ultraviolet divergences is not a very 
delicate issue here since we are dealing with a theory that has a natural 
cutoff. However, it is worth mentioning here how the removal of divergences 
can be implemented. For the initial conditions ( |3.24D , the divergences in 



< <p 2 (x, t) > can be shown to be the same ones obtained in the lowest order 
adiabatic approximation , i.e.: 

< <P 2 (x, r) > div = - [ [rf 3 k]-i— . (3.26) 

Introducing a cutoff at the physical momenta k v /r and k±, we can easily see 
that the above integral has both quadratic and logarithmic divergences. To 
renormalize the theory we absorb the quadratic divergences in the bare mass 
Xv 2 and the logarithmic ones in the coupling constant A. After a few simple 
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manipulations (that involve adding and subtracting the appropriate terms in 
(|3.10|) ) we can write the equation for \ as 

X(r) = *(7b) + ^($ 2 (r) - $ 2 (r )) + ^ / M 3 k]{|/ k (r)| 2 (l + 2 n k ) - 

- ^tt}-^/^]^^ ( 3 - 27 ) 

2^ k (r) - 1 r Z 7 w k (r ) 

where ct> k (r) = (k 2 /r 2 + k 2 i + x(to)) 1 / 2 and the renormalized coupling A r and 
Z are defined as: 

X = \ r /Z, Z = l-X r 5\, 5\ = — f[d 3 k]^--. (3.28) 

The initial value x( r o) comes from solving the gap equation ( |3.10|) at the 
initial time. As we discussed earlier it is useful to rewrite eq. ( |3.28|) in order 
to discuss triviality. Namely we have the inverse equation: 

A -db (3 - 29) 

The result is that at a fixed value of A there is a maximum renormalized 
coupling which decreases with the logarithm of the the cutoff. As we take 
the cutoff away, the renormalized coupling constant goes to zero, signifying 
the triviality of the theory as we remove the cutoff as first discussed by Baker 
and Kincaid ]23 and Bender et. al. |24 



To make our presentation more clear it may be useful to put together the 
set of equations we will solve. They are 

$i(r) + r-^r) + X (r) ^(r) = HS a , (3.30) 
X(r) = X(r ) + ^(r) - $ 2 (r )) + ^ J M 3 k]{|/ k (r)| 2 (l + 2 n k ) - 

-^}-^/M 3 k]^^ k , (3.31) 
2uj k {T)) r Z J LU k {r ) 

h + (^ + k 2 ± + X (r) + ^)/ k = 0. (3.32) 

It is worth noticing that in these equations the "bare parameters" A 
and v 2 do not appear. Thus, the equations are entirely written in terms 
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of renormalized quantities A r and xo- The value of xo is obtained from the 
gap equation (with the same cutoff in all the integrals). Once this is done, 
the value of the cutoff can be changed at will in ( |3.3U| )-( ]3.32| ) (provided we 
don't cross the Landau pole). 

For future use, it is convenient to compute the initial value of x, which is 
determined by equation (|3.31| ). It simply reads 

NX r r 3 n(k) 



X(r )(l + ^|[rf 3 k]-^) = 2A r $ i (r )<j. i (r ) 



2r J u^(t ) 



NX r r 3 n(k) k 



To J u k {r ) t -ujz 

(3.33) 

From here we can see that x( T o) has two contributions. The first term in 



the r.h.s. of ( |3.33| ) carries the contribution of the mean values: as expected, 
x( T o) is sensitive to the existence of initial velocities. The second term, which 
is always negative, arises from the finite temperature part of the initial state. 
It clearly shows how the expansion tends to reduce the initial value of x- 

It is instructive at this point to look at the simpler equations that arise 
in flat space for a spatially homogeneous expansion. In that case we have: 

$i(t)+x(t)$i(t)= H6 a , (3.34) 
x {t) = -Xv 2 + X^(t) + NX I [d 3 k] \ f k (t)\ 2 (l + 2 n k ), (3.35) 



f k + (k 2 + X (t))fk. (3.36) 
We also have the vacuum relations: 

ml = -Xv 2 + Xfl + NX[ [d\] 1 (3.37) 

J 2Jk 2 + ml 



Which yields: 

X (t) = ml + X$ 2 (t) - Xfl 

+NX [ [rf 3 k] [|/ fc (t)| 2 (l + 2 n k ) - 1 ]. (3.38) 
J 2^k 2 + ml 
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These equations govern the growth of instabilities. We see that there are un- 
stable modes fk whenever k 2 + x is negative so that only the long wave length 
modes can go unstable. If we ignore the quantum fluctuation contribution, 
X can go negative only when 



so that the most unstable case has $ = O.One also has that the bare coupling 
must satisfy 



for there to be any instabilities. Once instabilities grow, then they cause ex- 
ponential growth of the modes in the mode sum which contributes a positive 
quantity to the equation for \- Thus the stronger the coupling the quicker 
X returns to a positive value. In determining the parameters of our effective 
field theory we found that we are in the strong coupling regime (A r >> 1 
Thus we have the situation where any instabilities are quickly suppressed. 
This is not true at weak coupling such as found in the early universe prob- 
lems. In our approximation, we found that because of the triviality constraint 
the renormalized coupling constant is constrained to be less than or equal 
to about 10 is we use A = lGeV. We merely comment here that if we 
were able to consider higher values of A r then the relaxation of instabilities 
would be even faster and this would only dampen the possibility of producing 
disoriented chiral condensates. 

4 Initial Conditions, Correlations and Domains 
4.1 Reasonable initial conditions 

How can we study the possibility of forming large domains in which the pion 
field is disoriented? As we described above, the picture we have in mind 
is the following: After the collision, a "Baked Alaska" type configuration 
is formed. When viewed in the natural "boost invariant" coordinate system 
this configuration is described as a state characterized by homogeneous mean 
values $i(r) and by quantum and thermal fluctuations. The quantum state 



A/* < 0, 



(3.39) 




(3.40) 
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after the collision is clearly unknown and we will be forced to assume reason- 
able initial values for the parameters consistent with being near local thermal 
equilibrium in a disordered phase at a time r following the collision. Ideally 
we would like to be able to study the growth of inhomogeneous instabilities 
in real space and time. However to simplify our calculation we are assum- 
ing that the system evolves such that all expectation values just depend on 
the proper time. This rules out a detailed study of domain growth and we 
have to obtain information about the growth of domains indirectly from the 
quantum correlation functions. These are parametrized by the proper time 
evolving interpolating phase space number and pair densities. Thus, we will 
not attempt to describe here the process of formation of a "peculiar" quan- 
tum state (which is an interesting but very difficult issue). We will rather 
study the evolution of different reasonable initial states constrained to have 
expectation values which depend only on the proper time r, focusing on the 
possibility of the development of instabilities and on the potential growth of 
long range correlations. The question is how to choose the initial conditions 
defining a "peculiar" , but not entirely unrealistic initial state? Without any 
rigorous justification, we will assume an initial state that is a disoriented (or 
displaced) thermal state. Let us describe it in more detail. 

Let us forget for the moment about the expansion and consider a much 
simpler situation: thermal equilibrium at temperature T. This is character- 
ized by a value of x = Xt and by mean fields n = tyt = 0, a = o~t = H/xt- 
The value of xt, which is obtained by solving the gap equation given in 
the appendix, is positive and fixes the "effective mass" of the quasiparticles 
(at very high temperatures we have xt X r T 2 and a T oc H/T 2 , i.e. the 
explicitly broken symmetry is restored at high temperatures). 

In choosing an initial state to model the situation after the collision it 
is reasonable to assume the conditions of local thermal equilibrium above 
the critical temperature so that we start in the disordered phase where the 
chiral symmetry is unbroken. In particular, we can simply take the initial 
values o-(t ) = a?, 7r(r ) = 7Tt and x( r o) — Xt and turn on the expansion 
at T (the results of this simulation will be explained in detail later). One 
should also consider the possibility of exciting the initial state with some 
extra kinetic energy. For this we should consider an initial distribution of 
velocities &(tq) that kick the initial mean values making it change in time (it 
is worth noticing that the expansion destroys the initial equilibrium and the 
fields start to change even without any kick). Other initial conditions we will 
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consider are those in which the initial state is a "disoriented" equilibrium 
where ti 2 (t ) + cr 2 (r ) = o\ + tx\ but the state is pointing in the wrong 
direction (say, along tt°). In fact, the formation of such disoriented state 
is not at all unlikely: simple estimates show that the energy required to 
"tilt" the initial state orienting it along the 7r° direction (instead of the a 
direction) is only 10 MeV/fm 3 0. Another possibility is to change the 
initial value of the magnitude of the 0(4) vector making it different from its 
high temperature equilibrium one. For all the above cases, the initial value 
of x w iH always be determined by the solution of the initial gap equation (the 
same for its initial derivative) evaluated for an equilibrium configuration at 
a temperature T . 

However, we should point out that the initial value of x will always be 
restricted to be positive. As we said, xij) is the effective mass of the quasi- 
particles. Therefore, taking a negative initial value for x implies that we 
are "turning the effective potential upside down." In our view, this cannot 
be the consequence of forming a "peculiar" initial state but must rather be 
the consequence of the cooling mechanism, which is entirely produced by the 
expansion. In fact, starting with a negative xo is what is done when study- 
ing this problem by using the quench approximation: one starts with a hot 
initial state and let it evolve in the low temperature effective potential. It 
should be clear by now that this is drastically different from our approach. 
We will study the self-consistent evolution of x starting from a "hot" initial 
value and follow its evolution. 

4.2 Instabilities and correlations 

During the nonequilibrium evolution it is possible that, for certain time in- 
tervals, the value of x{ T ) becomes negative (we will discuss some examples 
below). When this happens there is an instability in the system. As clear 
from ( |3.32| ), if x < there are long wavelength modes that become unsta- 
ble: their amplitudes start growing exponentially (the factor l/4r 2 very soon 
becomes negligible). The existence of such unstable modes is the crucial in- 
gredient needed for the development of structure through the mechanism of 
spinodal decomposition. Let us briefly describe now how we will study this 
issue here. The intuitive picture is clear: the system cools as it expands 
evolving (in a fully nonequilibrium way) towards a stable low temperature 
state. In our case, such state (the vacuum) is characterized by the values 
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ff = 0, a — and x = m 1- The existence of an instability means that the 
homogeneous configuration is not energetically preferred and that any small 
inhomogeneity seed will tend to grow. As the 0(4) symmetry is explicitly 
broken, the growth of structure is only transient since in the long run the sta- 
ble state is again the homogeneous vacuum. The question one should try to 
answer is if during the nonequilibrium stage the instability is strong enough 
to form large domains in which the field is correlated and disoriented. If such 
domains do form, the correlations in the emitted pions can be detected. We 
should then examine the evolution of the correlation length (in the rapidity 
space) . 

A natural question that arises is how can one study the existence of long 
range correlations within our scheme, which is basically a mean field theory 
analysis. The textbook answer to this is that mean field theory can still be 
used to provide useful information about the typical scale of the correlations 



despite of the fact that the mean values are taken as homogeneous [2^]. In 
fact, this analysis enables us to compute the two point correlation functions 
which determine the behavior of the system when perturbed away from its 
homogeneous configuration. Thus, the correlation length obtained from the 
two point functions will characterize the growth of structure, at least in the 



linear regime. This point of view was used in [14, |15|). Our feeling however is 



that it is not so easy to interpret directly the two point correlation function. 
In the first order adiabatic basis we have: 



G°(x, y;r) = - f [d\y k ^-\- (1 + 2 n k (r) + 2Re[g k (r)e- i2 ^})(AA) 

We see that there are phases in the fourier transform G(k, r) between 
the number density n^{r) and the pair density g k {r) which can make the 
interpretation difficult. These interpolating operators only make physical 
sense when < 7r >= and for the stable modes. So at best we can look 
at the two quantities n k and g k and study their momentum dependence. 
If we have a particular model for these quantities which have proper time 
dependent masses as well as temperatures then can one extract correlation 
lengths (or mass scales) in a model dependent fashion. For example if the 
system expanded in local thermal equilibrium in a comoving frame we would 
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have: 



where /3(r) = 1/T(r) and £7 fc (r) = + x(r). 

The effective temperature could be calculated by first determining the 
hydrodynamical quantities e(r) and p{r) from the diagonal entries of the 
expectation value of the energy momentum tensor in the comoving frame as 



discussed in |29]. Then one determines the temperature (and entropy) from 
the relations : e + p = Ts ; de = Tds. 

We notice that this parametrization fails exactly when the pion mass 
goes negative which is the case we want to study. Also from our previous 
study of the production of pairs from strong electric fields in a related 1/N 
expansion [[^J , we know that the effective temperature in lowest order is not 
monotonically decreasing but is oscillating about a decreasing function with 
the plasma oscillation frequency. Thus until we add scattering (which occurs 
at next order in the 1/N expansion), the temperature parameter does not 
have this monotonicity property. At late times one expects after a period 
of entropy production that st = constant and T(r) = To (^) 17,3 as in the 
thermal equilibrium case. This is obviously one possible parametrization of 
the data in terms of two r dependent inverse length scales, the temperature 
and the mass of the pion. As we shall discover, this parametrization of the 
number density does not agree with our nonequilibrium evolution so that 
there are at least 3 proper time evolving length scales in this problem-the 
inverse mass of the pion, the inverse of the effective temperature and possible 
length scales describing domain growth. In what follows we will present our 
results for these two interpolating densities to see their general proper time 
development without making specific models in terms of various length scales 
except for the quasithermal model which does not reproduce the results. 

The computation of the pion two point function in general is straight for- 
ward. The two point function is obtained by inverting the inverse propagator 
matrix as discussed in the appendix. In general we get for the inverse of the 
pion two point function: 

G~l{x,y)ij = -S i:j p + x(x)]S(x -y)- iri(x)G 0xx (x - y)irj{y) (4.3) 

We see that only when tt^x) = can one obtain the Green's function directly 
in terms of the modes used to calculate the quantum field <fi. Also only in 
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that case do we get a direct interpretation of the Fourier transform of the 
Green's function in terms of the (proper) time dependent interpolating phase 
space number and pair densities described above. As discussed earlier it is 
simple to extract n and g from the solution / of the mode equation. 

Although the equal proper time correlation function depends only on the 
time dependent number and pair densities, the full non-equal time correlation 
function has more phase information and is given by (when < 7r >= 0): 

G°(x ± - x' ± ; rj - i/; r, r') = f [d^e'^'* '> (4.4) 

V r V T J 

{(1 + 2 n k )Re[r k (r)f k (r')} + 2Re\g k f k {r)f k {r')}} (4.5) 

where here n, g refer to their value at r . 

The only thing left to discuss is the choice of initial proper time r. One 
can estimate this quantity in a hydrodynamical model with Landau's choice 
of initial conditions (namely energy density of a Lorentz contracted pancake 
being given at an time zero). If one assumes that the the collision produces 
a quark-gluon plasma in equilibrium, one can run the hydrodynamical code 
for initial energy densities appropriate to the collision of Lorentz contracted 
disks and determine the proper time when the critical temperature is around 
200MeV or slightly above the chiral phase transition. We want to start in 
the quark - gluon plasma phase above the chiral phase transition which places 
constraints on the initial energy density present at r . We then assume that 
slightly above this transition it is reasonable to model the chiral transition 
with the effective Lagrangian of the sigma model. Not knowing exactly what 
this proper time is, we will here consider reasonable initial proper times 
(1/m/c < t < Afm/c and study the effect of the initial proper time r on 
the production of instabilities. Taking larger proper times as the starting 
point for our calculation reduces even further the possibility of instability 
growth so that our results present an upper limit on the growth of domain 
sizes in this model. 

5 Results 

We have performed numerical simulations on the connection machine CM-5 
using a grid that has lOOOOr modes. (We start at t = 1). The grid is 100 
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modes wide in the transverse direction and 100-^ modes in the rt direction. 

TO ' 

That is we choose: 



A dk n A 

ak± = , „ = . 

100 lfm 100 

The first issue we will examine here is the existence of instabilities. In pre- 
vious works the presence of unstable modes was assumed by imposing the 
quench approximation. In our numerical studies we find that, due to the 
strong coupling, that a wide class of initial conditions consistent with rea- 
sonable fluctuations found in a thermal distribution no instabilities develop. 
In our simulations the one thing we did not change was the value of the 
composite field x which was fixed to be the solution of the gap equation in 
the initial thermal state. We also maintained the constraint that the initial 
values of the expectation values of the field, namely vr(ro) and cr(r ) satisfied 
the relationship: 

7r 2 (r ) + a 2 (r ) = o 2 T 

where ot is the equilibrium value of $ at the initial temperature T. This last 
constraint is reasonable since it doe not cause much energy to make a rotation 
in the direction of 7T». Thus we chose different initial a and 7Tj consistent with 
the above constraints and then probed the effect of different initial values for 
a and 7Tj. Since the results for 7Tj 7^ were similar to those when a 7^ we 
mainly present here the results for the latter case. We have also surveyed 
other possibilities that violate the above constraints such as choosing initial 
a = 0. In that case we found no instabilities even with a 7^ 0. 

First let us consider the case when the initial value cx(to) is the thermal 
equilibrium one corresponding to a temperature of 200MeV. We varied the 
value of the initial proper time derivative of the sigma field expectation value 
and found that there is a narrow range of initial values that lead to the growth 
of instabilities. Namely 

.25 < \&\ < 1.3 

Surprisingly when \&\ > 1.3 instabilities no longer occur. 

Figures 1-2 summarize the results of the numerical simulation for the 
evolution of the system ( |3.30| ) (|3.32| ) . We display the auxiliary field x in units 



of fm~ 2 , the classical fields $ in units of jmT x and the proper time in units 
of fmT x (lfm" 1 = 197MeV). In Fig. 1 and Fig. 2 the proper time evolutions 
of the auxiliary field x field are presented for two different values of /„-, where 
the initial conditions were fixed at r = 1. The initial conditions in all the 
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Figure 1: Proper time evolution of the x field for four different initial condi- 
tions with f w = 92.5MeV. 

Figure 2: Same as Fig. 1, but for f n = l25MeV. 

cases are with the equilibrium value of the x held & t the corresponding initial 
temperature, where as the various initial conditions of the classical fields $ 
are chosen to satisfy 7? 2 (t ) + a 2 (r ) = a\ where <j\ is the equilibrium value 
of $ at the initial temperature T. In the case of f n = 92.5MeV we see that 
the instability lasts for less than 3/m. As discussed earlier, the regime of 
exponential growth of the unstable modes occurs whenever x < 0. We notice 
that when we rotate the expectation value from the a to the tt 1 direction 
initially then the instability regime has the same typical "survival" time. 
We also notice that if we choose the time derivative to be zero and start 
from an equilibrium configuration, then the expansion alone is insufficient 
to generate instabilities. We find that in order to generate instabilities we 
require fluctuations in the classical kinetic terms such as & or -hi and as 
discussed earlier these initial conditions must be in a very narrow range 
to produce instabilities. Thus the rapid quench conditions assumed by other 
authors comprises only a small region of the phase space of initial fluctuations 
expected in an initial thermal distribution. As expected, at late proper times, 
the auxiliary field approaches its equilibrium value of m 2 . For the case when 
f n = 125MeV, a value favored by fitting the low energy scattering, we see 
a very similar behavior, showing that our main results concerning the small 
range of initial conditions that lead to instabilities are not affected by a 30% 
change in the value of one of our parameters. 

A crude estimate for the size of a disoriented chiral condensate is to 
multiply a typical survival time length for the instabilities by the speed of 
light. If this estimate is valid then the size of these regions are of the order 
of a few fermi. Of course one needs to study the growth of inhomogeneous 
instabilities to make any definite statements about this size. We were hoping 
that the correlation functions we have calculated could be interpreted easily 
in terms of a length parameter associated with the size of these regions. 
However as we will see below such a naive hope is not to be satisfied. In 
Fig. 3 we plot the proper time evolution of the classical fields a and n for 
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Figure 3: Proper time evolution of the a and ir fields for the following initial 
conditions: Solid line is for er(l) = a T , tt 1 (1) = and cx(l) = 1. Dotted 
line is for er(l) = a T , tt 1 (1) = and <j(l) = 0. Dashed- dotted line is for 
(7(1) = a T , tt^I) = and = -1. Dashed line is for cr(l) = 0, 7r*(l) = a T 
and d(l) = -1. At T = 200Me1/, a T = O.Sfm- 1 for / w = 92.5MeV and 
a T = 0.5/m- 1 for / w = 92.5Me\/ 

Figure 4: Proper time evolution of the \ field for three different initial proper 
times r = 1,2.5,4 with f w = 92.5MeV, a(T in ) = a T , ir l (T in ) = and 
cr(r in ) = -1. 

the same two choices for /„. as in Fig. 1 and Fig. 2. We see that the sigma 
field asymptotes to its vacuum value /„- and the 7r field gradually converges 
to its equilibrium value of zero. In Fig. 4 we show the effect of starting 
the initial value problem at later times, namely r = 2.5, 4/m and compare 
them to the case previously studied with r = lfm which had a modest 
region of instability growth. We see that as we increase To we decrease the 
possibility of instability growth and by these late times it is not possible 
to produce instabilities even with kinetic energy fluctuations. In Fig. 5 we 
study the effect of the initial temperature on our time evolution problem. 
For f w = 92.5 the critical temperature is 160MeV in the absence of explicit 
symmetry breaking. We see that in the vicinity of the critical temperature 
the effects of varying the initial temperature is minor. In Fig. 6 we study 
the proper time evolution of the effective number density for various initial 
conditions. In thermal equilibrium one expects for an isentropic expansion in 
boost invariant coordinates that sr = constant. Since the number density is 
proportional to the entropy density one expects that once particle production 
stops that n(r)r — > constant. We see this trend in this figure, showing that 
we are reaching the out regime as the system expands. The breaks in some of 



Figure 5: Proper time evolution of the x fi e ld for three different initial 
thermal distributions with T = 200, 164, 150MeV for the initial conditions 
v(r in ) = o- r , -K l {n n ) = and &(r in ) = -1. 
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Figure 6: Proper time evolution of the produced particle density nr = 
dN/dr]dx± for the same evolution shown in Fig. 1 

Figure 7: Slices of k v = and p = |kj_| = of the proper time evolution of the 
interpolating phase space particle number density n(k v , kj_, r) for er(l) = cr T , 
7r J (l) = and <t(1) = —1 compared with the corresponding local thermal 
equilibrium densities riT(k v , kj_, r). 

these curves at early values of r are a result of the fact that the interpolating 
number density cannot be defined when x is negative. Next we are interested 
in knowing how our results differ from the case where the system evolves in 
local thermal equilibrium which is described by two correlation lengths, the 
inverse of the effective pion mass associated with x > an d the inverse of 
the proper time evolving effective temperature T(r) = To( — ) 1 ^ 3 discussed 
earlier. We see from Fig. 7 that in the case that er(l) = <jt, 7t'(1) = and 
<t(1) = — 1, where maximum instability exists, complex structures are formed 
as contrasted to the local thermal equilibrium evolution. The interpolating 
phase space distribution n(k v , kj_, r) obtained numerically, clearly exhibits a 
larger correlation length in the transverse direction than the equilibrium one 
and has correlation in rapidity of the order of 1-2 units of rapidity. We notice 
that in both directions there is structure which does not lend itself to a simple 
interpretation. On the other hand the local thermal equilibrium evolution 
is quite regular apart from the normalization of the distributions that are 
changing with time due to oscillation in the quantity x( T ) which is damped 
to its equilibrium value once the system expands sufficiently Other authors 
have suggested that it is possible to extract the sized of the domains of DCC's 
from the coordinate space correlation function G(x,y,r) by considering the 
width of the distribution as a measure of the size of the domain. However, it 
is clear from our detailed numerical results that there are several length scales 
affecting the momentum space distribution apart from the effective pion mass 
and temperature. The spatial Green's function depends on not only n k and 
gk but also the phases y& and would be even harder to interpret than n and 
g separately. In Fig. 8 we look at the time evolution of the interpolating 
number operator for the case when we did not have any fluctuations in the 
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Figure 8: As in Fig. 7, but for tr(l) = a T , 7r*(l) = and <r(l) = 0. 

Figure 9: Slices of A;^ = and p = |kj_| = of the proper time evolution 
of the real part of the phase space pair density g(k v , k^, r) for cr(l) = cry, 
7r*(l) = and <r(l) = 0, -1. 

kinetic energy. Here we are still far from equilibrium and we see that unlike 
the equilibrium case the correlation length in rapidity space is not decreasing 
with proper time. In this case where there are no instabilities we do not 
see complicated structures. The transverse distribution is similar to the 
equilibrium case. The pair density function g(k v , k±, r) is even more elusive 
to parametrize than the single particle distribution function. In the case 
where there are no instabilities (<j(l) = 0) we see in Fig. 9 that although the 
transverse distribution is relatively simple, the distribution in the 77 direction 
(whose fourier transform gives the rapidity distribution) has many length 
scales. When we have instabilities (<r(l) = —1) then both distributions are 
complicated and possess several length scales as seen in the lower part of 
Fig. 9. 

6 Conclusions 

In this paper we have performed numerical simulations in the regime of the 
chiral phase transition in the linear sigma model for a wide variety of initial 
conditions starting above the critical temperature for an expanding plasma of 
pions and sigmas. Assuming that this model gives a reasonable description in 
this temperature range, we found initial conditions where instabilities grew 
in the scenario required for the formation of disoriented chiral condensates. 
The constraints on our model which are imposed in order to approximate the 
low energy physics of pion interactions however require a large renormalized 
coupling constant. This had the effect of rapidly damping the instabilities 
and we found no evidence in this model for large domains of disoriented 
chiral condensates. We did, however see rather large departures in the phase 
space number density from one which would result from an evolution in local 
thermal equilibrium. These departures show a narrowing of the momentum 
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space distribution which could be interpreted as a larger spatial correlation 
length. However we found no simple method of extracting correlation lengths 
from our results. In our simulations we assumed a reasonable mechanism for 
cooling- namely the expansion of the plasma following its production in a 
collision. However we did ignore scattering, which occurs at next order in the 
1/N expansion. We also did not allow for general inhomogeneous fluctuations 
which would allow us to actually study the growth of disoriented domains. 
This would require a study of the growth of inhomogeneous perturbations 
about our homogeneous background and would involve a much more difficult 
numerical computation. We hope to carry out such a computation in the 
future. It is also possible that another model, such as the non-linear sigma 
model might better describe the dynamics in this regime and might lead to 
different conclusions from those found here. In the future we also hope to 
include scattering effects which will introduce another time scale, namely 
the equilibration time scale into the problem, as well as study the nonlinear 
sigma model to see if the results differ significantly from those found here. 



7 Appendix I- a model in the large N ap- 
proximation 

The 0(4) a model is described by the Lagrangian 

L = {-d& ■ <9$ - -A($ • $ - v 2 ) 2 + Ha}. (7.1) 

where the mesons are organized in an 0(4) vector $ = (jr, a) . The counting 
of the large N expansion is made explicit []T7| by introducing a composite 
field x defined by % = A($ • $ — t> 2 ). It is easy to show by appropriate 
rescalings that the large N expansion is obtained by integrating out the $ 
field and then performing a steepest descent calculation on the remaining 
X path integral. That is the large N expansion is a loop expansion in the 



composite field x propagator [19]. Thus we add to the above Lagrangian the 
constraint equation: 

[ x -A($-$-^)] 2 

4A ' 1 ] 
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which exactly cancels the quadratic term. We thus obtain the alternate 
Lagrangian: 

L 2 = -^(D + + ^ + l - X v 2 + Ha. (7.3) 

The generating functional of the Green's functions is given by the Path 
Integral: 

Z[J,S] = j D X D<$> exp {i Jd 4 x[L 2 + J^ + Sx]}] = exp[ i W(J, X )}.(7A) 

We can now perform the Gaussian path integral over the $ field. Evaluat- 
ing the remaining x integral at the stationary phase point of the resulting 
effective action, and then Legendre transforming: 

x] = W[J, S]-J d 4 x[J(x) ■ $(x) + S(x) X (x)}, (7.5) 

we obtain the lowest order (in large N) result: 

r[$, X ] = J d A x[L 2 ^,x,H) + ^iVtr InGo 1 ], (7.6) 

where 

G^(x,y)=t[D + x(x)}6\x-y). (7.7) 

For the 0(4) sigma model N = A. From this effective action we immediately 
get the equations of motion for the fields $ and the equation of constraint 
(gap equation) for the composite field x '■ 

p + x (x)]m = p + x (x)]a = H, (7.8) 

x = -Xv 2 + X(a 2 + 7T ■ tt) + XNG (x, x). (7.9) 

In the static case when we are considering symmetry breaking, we obtain 
lowest energy state when: 

X a = H; x = ~^v 2 + Xa 2 + XNG . (7.10) 
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By considering the fields ($, x) to be a 5 dimensional vector and the 
matrix inverse propagator to be: 

G- 1 = (7.11) 

we obtain by inverting this matrix the relevant it, a, and \ propagators. 
Performing this inversion we obtain for the vacuum Feynman x inverse prop- 
agator: 

GUP 2 ) = Go^P 2 ) + (7.12) 

where 

Go \ x {p 2 ) = ^ + fn(p 2 ) (7.13) 

is the inverse propagator in the absence of symmetry breaking and the po- 
larization II = iGl is given by 

n( P 2 ) = -iJ [A] (x - q 2 Y\x - (p + qfr 1 - (7.14) 

In our simulations we have a three dimensional cutoff. If we integrate 
over q° and perform the remaining three dimensional integral with a cutoff 
A we obtain explicitly: 



rV)= ! 



8tt 2 



ln(- + Wl + -) + -,! ---In 



1 / T 1 / W Ml- ^)(1+^ 



v 



1 V v 2 ' 2\ /r \ l • (1 . . 1^-1(1 • .r 



(7- 



where x = m/A and m is the pion mass: x — m2 - 

For the pion inverse propagator we obtain in the vacuum sector: 



G-liP% = ^{P 2 ~X) = UP 2 - m 2 ). (7.16) 

Thus we see that x — m2 is the pion mass squared. In our initial value 
problem one can have 7Tj ^ 0. For that case (or when there are external 
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isospin sources present), the pion inverse propagator becomes for constant 
external sources: 

G^Vk = 5 l3 (p 2 -x)- vr.Go xx (P>r (7-17) 
For the a inverse propagator we obtain: 

G„ 1 Ap 2 )=p 2 -m 2 -a 2 G 0xx (p 2 ). (7.18) 
Thus the a mass m 2 is determined by the relationship: 

ml = m 2 + a 2 Re[G 0xx (ml)}. (7.19) 
The axial current in this model is given by : 

^(x) = [7i l d,a - ad,n% (7.20) 
The PCAC condition is given by: 

d.A^x) = Hn\x). (7.21) 

Therefore one has 

H = Urn 2 - (7.22) 

Since we also have that the vacuum is defined by x a — m2(J — H, we 
immediately find that a — f n . We can rewrite the gap equation in terms of 
physical quantities as: 

m 2 = -\ v 2 + Xf 2 + AiVGo (x, x) , (7.23) 

where 

This relationship allows us to determine the bare mass — Xv 2 in terms of the 
renormalized mass m, the pion decay constant f w and the cutoff A 

To determine A r we must look at either the a mass which is not well 
determined by experiment, or the low energy scattering data such as the s 
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wave scattering 1=0 scattering amplitude. Although at tree level the sigma 
mass is an arbitrary parameter, that is not true for the quantum theory. The 
reason for this is that we only want to consider this model as an effective field 
theory with a cutoff, so that the bare coupling constant is always positive. 
(Renormalized ordinary perturbation theory requires that the bare coupling 
constant is negative as we take away the cutoff). Once there is a cutoff, then 
the property of the exact 0(N) sigma model, as determined by lattice or 
strong coupling expansions is that the renormalized coupling is a monotoni- 
cally increasing function of the bare coupling reaching a finite maximum as 
the bare coupling goes to infinity. In four dimensions, this maximum value 
decreases logarithmically with the cutoff. This exact feature of the O(N) 
sigma model is preserved in the large N approximation which is a virtue of 
our approximation. However for reasonable cutoff (say 1 Gev) this puts and 
upper bound on the a mass of around 3m as we shall see below. The low 
energy scattering data presents another problem. This data is quite difficult 
to fit in perturbation theory because perturbative results violate partial wave 



unitarity. Basdevant and Lee |22j were able to fit the s wave data by first 
calculating to one loop and then enforcing unitarity by doing a Pade approx- 
imant. Our 1/N approximation is related to the Pade approximant since it 
sums the bubbles of the perturbative result, but it does not automatically 
obey partial wave unitarity. Nevertheless it is important to see what values of 
A r give reasonable s wave phase scattering amplitudes in our approximation. 
Let us therefore turn our attention to determining the low energy scattering 
amplitude in the large N limit. The pion-pion scattering invariant T-matrix 
is given by: 

T = dijdkiAis) + 5 ik 8jiA{t) + 8 u 5 jk A(u) (7.25) 

where the isospin indices i, j, k, I are coupled to the four momenta pi,P2, Ps, Pi 
such that s = (pi + P2) 2 , t — {p\ — p-s) 2 u — (pi — Pa) 2 ■ In leading order in 
large N the amplitudes A(s), A(t) and A(u) are exactly the x propagator. 
Namely 

A(s) = -G xx (p 2 = s). (7.26) 

The large N expansion preserves the current algebra so that the usual low 
energy theorem is exact. That is, for small s,m 2 one easily shows that 

A(s) 

■1 7T 
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This amplitude is independent of the coupling constant A. We thus need to 
go above threshold to determine the coupling constant. The 1 = scattering 
amplitude is 

A = 3A(s) + A{t) + A(u). (7.27) 

The s-wave scattering amplitude is obtained by integrating the the 1=0 scat- 
tering amplitude over angles. 



i / Am 2 r l 
f l=0 = e*Msin 6(s) = -^Jl -— dz A , (7.28) 

62tt V s J -i 

where z = cos 9 and 9 is the scattering angle in the s channel center of mass 
system. 

It is useful to describe the theory in terms of the renormalized coupling 
constant A r which depends on A as well as the cutoff A. The running renor- 
malized coupling constant is determined by the renormalization group invari- 
ant composite field propagator G x x .We choose to define the renormalized 
coupling constant A r as the running coupling constant at q 2 = 0, for the 
unbroken mode of the theory. That is: 

Ar = iTliH (7 ' 29) 

Rewriting the propagator in terms of A r ,one immediately gets a finite expres- 
sion for the running coupling constant. 

o \ 

2A r (g 2 ) = G x x (q 2 ) = ^-^ , (7.30) 

1 + VVIW) - 

with 

n r (g 2 ) = n(g 2 )-n(g 2 = 0). (7.31) 

In Fig. 11 we plot the real part of the s wave scattering amplitude for 
A r = 7.3, = 125MeV and A r = 7.3, f n = 92.5MeV. By comparing these 
curves to the unitarized perturbation theory results of Basdevant and Lee 



[2^| we find reasonable agreement in the regime 2 < ys/m < 2.6 for f w = 
125MeV .This value of /„- is obviously the preferred one if we want a closer 
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Figure 10: The real part of the partial wave amplitude fi=o for 1 = com- 
puted from the linear cr-model in the lowest order of the 1/N expansion for 
/„- = 92.5MeV and for /„. = 125MeV. The value of A r used was 7.3. 



agreement of our results with those of Basdevant and Lee. (However since 
in our large N expansion we have only a modest value of N, namely N = 4 
we might expect 25% corrections at next order in 1/N). We notice the 
existence of the sigma resonance at a mass of around 3m. These values of 
f n and m a are in reasonable agreement with the fits of Basdevant and Lee 
using their Pade analysis of ordinary perturbation theory. We see that in 
our approximation there is a slight break down in s-wave unitarity near the 
peak, (s- wave unitarity requires that the real part of this amplitude to be less 
than 1/2). This breakdown in s-wave unitarity occurs in the leading order 
in large-N approximation for renormalized couplings which are larger than 
three. 

One can approximately calculate the sigma mass from the zero of the real 
part of the inverse propagator. In terms of the renormalized coupling defined 
above we have: 

9X f 2 

< - '" 2 + ft| i+ww'' (7 ' 32) 

This equation, which leads to a similar value for the mass of the sigma as 
found from the peak in the scattering, shows that there is an upper bound 
on the sigma mass which depends on the chosen cutoff, since A r decreases 
monotonically with the cutoff. 

The phase structure of the a model in this approximation has been studied 
extensively by W. Bardeen and Moshe Moshe | 20| . For zero temperature the 



expectation value < T 00 > in our initial density matrix is the correct energy 
density functional to study. At finite temperature it is the free energy density 
expectation value that is needed to determine the phase structure. Explicitly 
at zero temperature one has that the energy density functional is: 

M pV2A2 \ AT „A2 

W(* ■ S, x) = 7rr- 2 X 2 H- ] + j [* ■ $ - « 2 - T7- 2 XH—)} 2 - (7-33) 
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Minimizing this equation with respect to x gives the gap equation 

x = A(*.*-£)-^ln(^), (7.34) 

which then gives x as a functional of $ • $. This then implicitly determines 
the energy functional as a functional of only Bardeen and Moshe Moshe 
point out that if one uses the endpoint solution x — below the minimum 
of the potential, then one gets a real "effective potential". 

At finite temperature the phase structure of the cutoff theory is deter- 
mined instead by the free energy density functional: 

rx,/* t \ iV 2 .e^A 2 . AT [x/t 2 r dF(y), 



X 48 Jo dy 

+ —T 2 F(- 

24 v : 

(7.35) 



^ - 9 AT , 9 AT , ,eA 2 A^ 9 . Y . l9 

+ - $ • $ - v 2 + -A -v n + — T F 4 

4 L 16tt 2 16tt 2A x 24 V T 2 ^ J 



where 



The gap equation is now: 

X = -^l»(f> + ^<£). (7.37) 

Using this free energy one finds that the critical temperature is determined 
from: 

NT 2 

^ = & ( 7 - 38 ) 
These relations are used to start our calculation above T c . 
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